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Abstract. We introduce a generalized Ulam method and apply it to symplectic dynamical maps with a 
divided phase space. Our extensive numerical studies based on the Arnoldi method show that the Ulam 
approximant of the Perron- Frobenius operator on a chaotic component converges to a continuous limit. 
Typically, in this regime the spectrum of relaxation modes is characterized by a power law decay for small 
relaxation rates. Our numerical data show that the exponent of this decay is approximately equal to the 
exponent of Poincare recurrences in such systems. The eigenmodes show links with trajectories sticking 
around stability islands. 

PACS. 05. 45. Ac Low-dimensional chaos - 05.45.Pq Numerical simulations of chaotic systems - 05.45.Fb 
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1 Introduction 

The properties of two-dimensional (2D) symplectic maps 
with dynamical chaos have been studied in great detail 
during last decades both on mathematical (see e.g. [TJ[2] 
and Refs. therein) and physical (see e.g. [3,4,5 and Refs. 
therein) grounds. A generic and nontrivial behavior ap- 
pears in maps with divided phase space where islands of 
stability are surrounded by chaotic components. A typical 
example of such a map is the Chirikov standard map [3j[4] 
which often gives a local description of dynamical chaos 
in other dynamical maps and describes a variety of physi- 
cal systems (see e.g. [6]). This map is characterized by one 
dimensionless chaos parameter K and two dynamical vari- 
ables y which have a meaning of phase and conjugated 
action: 

K 

y = y + — sin(27rx) , x = x + y (mod 1) . (1) 

Here bars mark the variables after one map iteration and 
we consider the dynamics to be periodic on a torus so that 
< x < 1, < 2/ < 1. 

For small values of K the phase space is covered by in- 
variant Kolmogorov-Arnold-Moser (KAM) curves which 
restrict dynamics in action variable y. With the increase 
of K more and more of these KAM curves start to be de- 
stroyed and above a certain K c all curves disappear and 
dynamics in y becomes unbounded. In 1979 Greene [7 
argued that the last KAM curve has the golden rotation 
number r = r g = ((x t — xo)/t) = (\/5 — 1)/2 with the criti- 
cal Kg = 0.9716... (here t is given in number of map itera- 
tions; there is also symmetric critical curve at r = 1 — r g at 
Kg). A renormalization technique developed by MacKay 



[8] allowed to determine K g = 0.971635406 with enormous 
precision. The properties of the critical golden curve on 
small scales are universal for all critical curves with the 
golden tail of the continuous fraction expansion of r for 
all smooth 2D symplectic maps H . Further mathematical 
[9] and numerical [10] results showed that the actual value 
of K c is indeed very close to K g (K c — K g < 2.5 x 10 -4 
according to [10]) and it is most probable that K c = K g . 

The results of Greene and MacKay [7,8 gave a funda- 
mental understanding of the local structure properties of 
symplectic maps in a vicinity of critical invariant curves 
but the global properties of dynamics on a chaotic com- 
ponent still keep their mysteries. For K > K g the golden 
KAM curve is replaced by a cantori [TT] which can signif- 
icantly affect the diffusive transport through the chaotic 
part of the phase space [T2|fT3] . In addition there are other 
internal boundaries of the chaotic component with critical 
invariant curves which can affect statistical properties of 
chaotic dynamics. One of such important properties is the 
statistics of Poincare recurrences P(t) which is character- 
ized by a slow algebraic decay in time being in contrast 
to an exponential decay in a homogeneously fully chaotic 
maps (see [14,15,16,17,18,19,20 and Refs. therein). This 
algebraic decay P(t) oc \jt^ has (3 ~ 1.5. Such a slow 
decay appears due to trajectory sticking near stability is- 
lands and critical invariant curves and leads to even slower 
correlation decay with a divergence of certain second mo- 
ments. A detailed understanding of this phenomenon is 
related to global properties of dynamical chaos in 2D sym- 
plectic maps and is still missing. 

With the aim to analyze the global properties of chaotic 
dynamics we use the Ulam method proposed in 1960 [2T] , 
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In the original version of this method the phase space 
is divided in Nd = M x M cells and n c trajectories are 
propagated on one map iteration from each cell j. Then 
the matrix Sij is defined by the relation Sij = 7iij/n c 
where riij is the number of trajectories arrived from a 
cell j to a cell i. By the construction J2 i S^ = 1 and 
hence the matrix Sij belongs to the class of the Perron- 
Frobenius operators (see e.g. [22]) and can be considered 
as a discrete Ulam approximate of the Perron-Frobenius 
operator (UPFO) of the continuous dynamics. According 
to the Ulam conjecture [21] the UPFO converges to the 
continuous limit at large M. Indeed, this conjecture was 
proven for one-dimensional (ID) homogeneously chaotic 
maps [23 . Various properties of the UPFO for ID maps 
have been studied in [24,25 ,26 and further mathematical 
results have been reported in [27,28,29,30 with exten- 
sions to 2D maps. It was also shown that the UPFO can 
find useful applications in studies of dynamics of molecu- 
lar systems [31] and coherent structures in dynamical flows 
[32] . Recent studies [33,34 traced similarities between the 
UPFO, the corresponding to them Ulam networks and the 
properties of the Google matrix of the world wide web net- 
works. 

While for homogeneously chaotic systems the Ulam 
method is well convergent to a continuous limit it is also 
well known that in certain cases the discretization leads to 
violent modifications of system properties (see e.g. [28]). 
For example, for 2D maps with a divided phase space the 
UPFO destroys all KAM curves and thus absolutely mod- 
ifies the system properties (see e.g. discussion in [33 ). The 
physical origin of these unacceptable modifications is re- 
lated to a small noise, introduced by the coarse-graining, 
which amplitude is proportional to the cell size 1/M. This 
noise allows trajectories to penetrate through invariant 
curves leading to a broadly known opinion that the Ulam 
method is not applicable to the Hamiltonian systems with 
divided phase space. 

In this work we show that the Ulam method can be 
generalized in such a way that it becomes applicable to 2D 
symplectic maps with a divided phase space. We use this 
generalized Ulam method to investigation of the Chirikov 
standard map at the critical parameter K g and at large 
values of K when the phase space has small stability is- 
lands. Our extensive numerical simulations allow to obtain 
new features of the global chaotic dynamics in such cases. 
We also show that this method can be applied to other 
maps, e.g. the separatrix map or whisker map [4 . 

The paper is constructed as follows: in Section 2 we 
describe the generalized Ulam method and demonstrate 
its convergence for the map (fll) at K = K g , in Section 
3 we describe the Arnoldi metnod which allows to study 
the spectral properties of the UPFO in the limit of large 
matrix size up to Nd ~ 10 6 . The spectral properties of 
the UPFO are analyzed in Section 4 for the map (fil) at 
K = Kg and in Section 5 at K = 7. The case of the 
separatrix map with the critical golden curve is studied 
in Section 6, the discussion of the results is presented in 
Section 7. 
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Fig. 1. (Color online) Density plots of the eigenvector ipo 
of the UPFO with eigenvalue A = 1. The UPFO is ob- 
tained by the generalized Ulam method with a single trajec- 
tory of 10 12 iterations of the Chirikov standard map |l]) at 
K — Kg — 0.971635406. The phase space is shown in the area 
< x < 1, < y < 1/2; the UPFO is obtained from M x M/2 
cells placed in this area. The value of M for the panels is 25 
(first/top left), 50 (first/top right), 140 (second/bottom left), 
280 (second/bottom right), the corresponding dimension of the 
UPFO matrix S is N d = 177, 641, 4417, 16609 respectively. The 
probability density of the eigenstate is shown by color with 
red/grey for maximum and blue/black for zero. 

To make the Ulam method to be applicable for the 
symplectic maps with divided phase space we use the fol- 
lowing generalization of the method which we explain on 
an example of the Chirikov standard map Q. The whole 
phase space 0<x<l,0<y<lis divided on M x M 
equal cells. One trajectory is taken in the chaotic compo- 
nent (e.g. at xq = 0.1/27T, y$ = 0.1/27r) and is iterated on 
a large number of map iterations t, e.g. t = 10 12 . Then 
the UPFO matrix is defined as Sij = nij/J2i n ij where 
riij is the number of transitions of the trajectory from a 
cell j to a cell i. By the construction we have Sij = 1 
and hence this UPFO Sij belongs to the class of Perron- 
Frobenius operators. In this construction a trajectory vis- 
its only those cells which belong to one connected chaotic 
component. Therefore the noise induced by the discretiza- 
tion of the phase space does not lead to a destruction of 
invariant curves, in contrast to the original Ulam method 
[2T] which uses all cells in the available phase space. Since 
the trajectory is generated by a continuous map it cannot 
penetrate inside the stability islands and on the physical 
level of rigor one can expect that, due to ergodicity of dy- 
namics on one connected chaotic component, the UPFO 
constructed in such a way should converge to the Perron- 
Frobenius operator of the continuous map on a given sub- 
space of chaotic component. 

A mathematical prove of such a generalized Ulam con- 
jecture of the convergence of the UPFO built from one 
trajectory is not an easy task. Therefore, we performed 
extensive numerical simulations which confirm the con- 
jecture. With this aim we checked that the results for the 
spectrum and eigenstates of S remain stable when t is 
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changed from t = 10 10 to 10 12 , when we take another tra- 
jectory in the same chaotic component, and when the size 
M is increased (see detailed discussion below). To reduce 
the matrix size of Sij we use the symmetry property of 
the map ([!]) which remains invariant under the transfor- 
mation £— ^1 — x, y —> 1 — y so that we can consider cells 
only in the lower half square with 0<x<l,0<y<l/2 
which contains M 2 /2 cells. At K = K g we find that the 
number of cells visited by trajectory in this half square 
scales as N d « C d M 2 /2 with C d « 0.42. This means that 
the chaotic component contains about 40% of the total 
area that is in a good agreement with the know result of 

a- 

We used values of M in the range 25 < M < 1600. 
To be more precise, for practical reasons, we determined 
the UPFO (actually the integer numbers mj) for the two 
largest values M = 1600 and M = 1120 by iterating 
a single trajectory as described above and for smaller 
values of M we used an exact renormalization scheme 
by merging four neighbored cells (for a certain value of 
M) into one single cell (for M/2). In this way we ob- 
tained in an efficient way the UPFO also for smaller values 
M = 800, 560, . . . , 35, 25 without the necessity to reiter- 
ate the same classical trajectory. 

For t = 10 12 and M = 1600 we have about n c « 
2t/(C d M 2 ) w 1.8 x 10 6 transitions for each cell. This num- 
ber is rather large and relative statistical fluctuations are 
on a small level of \j yfn~ c ~ 10 -3 . 

The direct exact diagonalization of the matrix S can be 
done by standard computer routines which require mem- 
ory resources of iVj ~ M 4 double precision registers. The 
computational time scales at 7Vj ~ M 6 . Thus, for the 
map at K = K g we are practically limited to M = 280 
(with N d = 16609) as the maximum size for the full diag- 
onalization. At such M the statistical error is of the level 
rsj 10 -4 . Larger values of M can be reached by the 
Arnoldi method as it is discussed in the next Section. 

The eigenvalues Xj and corresponding right eigenvec- 
tors ipj(i) are defined from the equation 

N d -1 

^ Smi^jij) = Xj^j(m) . (2) 

2=0 

According to the Perron- Frobenium theorem [22] we have 
the maximal eigenvalue Ao = 1 with the corresponding 
eigenstate il>o(i) shown in Fig. [I] for four values of M. All 
values ^o(i) are non- negative in the agreement with the 
theorem and have the meaning of the probabilities in a 
given cell i. With the increase of M the state ipo(i) con- 
verges to a homogeneous ergodic measure on the chaotic 
component. The stability islands are well incorporated in- 
side the chaotic component. 

Another confirmation of the convergence of the UFPO 
in the limit of large M is presented in Fig. [2] In a first 
approximation the spectrum A of S is more or less ho- 
mogeneously distributed in the polar angle cp denned as 
Xj = \Xj\exp(i(pj) (see left column of Fig. [2]). The two- 
dimensional density of states p(X) clearly converges to a 
limiting curve. This density of states is normalized by 



J p(X)d 2 X = 1 (for a full spectrum of N d eigenvalues). 
It drops when |A| approaches to 1 but even at |A| « 0.9 
the convergence to a limiting curve is clearly seen. This 
is also confirmed by data with 400 < M < 1600 obtained 
from the Arnoldi method (which corresponds to a partial 
spectrum of 3000 — 5000 <C N d eigenvalues with largest 
| Xj | and is therefore not properly normalized) . 

The convergence of p(X) at Nd —> oo implies that the 
spectrum has a usual dimension d/2 = 1 corresponding 
to the dimension of the phase space. We note that the 
situation becomes different for dissipative maps where the 
fractal Weyl law determines the number of states in a 
given area of A that grows slower than N d (see [33| f35] 
and Refs. therein). Our direct computation of the number 
of states N\ in an interval 0.1 < A < 1 gives a linear 
dependence N\ oc N d . 

The properties of Aj, with |A| being close to 1 (see e.g. 
right top panel of Fig. [2j, and their scaling with M will 
be discussed in next Sections after a description of the 
Arnoldi method which is especially efficient in the compu- 
tation of such eigenvalues. 

Let us note that our special checks show that the vari- 
ations of Xj with the change of initial trajectory or its 
length t remain on the level of statistical accuracy 1 / ^Jn~ c . 
In the following we present data obtained with the tra- 
jectory length t = 10 12 which is close to a maximal com- 
putational effort used in [18 where t < 10 13 was used 
for the computation of the Poincare recurrences. Such a 
large value of t allows for the trajectory to penetrate into 
the very close vicinity of the critical invariant curve that 
becomes important at large M. 

3 Arnoldi method 

In order to capture features of small phase space structures 
(such as small stable islands) and to get a better approx- 
imation of the continuous limit (M —> oo) it is of course 
desirable to increase M further than the value M = 280 
accessible by the exact diagonalization. Fortunately the 
matrix S is very sparse with an average number of non- 
zero connecting elements per row (or per column) being 
k c « 5 (and maximum number of links ft m = 6, at K g ) and 
12 (and n m = 20 at K = 7). The value of a maximal num- 
ber of non-zero elements is determined by a local stretch- 
ing given by the monodromy matrix, thus we have ap- 
proximately k c ~ exp(/i) where h is the Kolmogorov-Sinai 
entropy (see [4] and discussion in [33]). Since k c <C N d we 
can calculate and store the matrix S for larger values of 
M and also effectively compute the product of S with an 
arbitrary vector with k c x N d operations. 

Furthermore, we are primary interested in the part 
of the spectrum with eigenvalues of modulus \Xj\ close 
to 1, or in other words with minimal decay rates jj = 
— 21n(|Aj|), in order to capture the long time properties 
of chaotic dynamics with the UPFO iterations. 

We have therefore used the Arnoldi method [36] which 
is perfectly adapted for this situation. This method is sim- 
ilar in spirit to the Lanzcos method, but is adapted for 
non-hermit ian or non-symmetric matrices. It has allowed 
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Fig. 2. (Color online) Spectrum Xj of the UPFO of the map 
at K = K g . First row : The left panel shows the eigenvalue 
spectrum in the complex plane for M — 280 and Nd — 16609 
by red/grey dots. The small blue/black square close to the 
region A = 1 is shown in more detail in the right panel with 
eigenvalues as red crosses. The green/grey curve represents the 
circle |A| = 1. Second row : In the left panel the Ritz eigen- 
values (blue/black squares), obtained by the Arnoldi method 
for M — 280 and with the Arnoldi dimension ua — 1500, 
are compared with the exact eigenvalues (red/grey dots). The 
right panel shows the modulus of the differences between the 
exact eigenvalues and the Ritz eigenvalues as a function of the 
level number j with eigenvalues sorted by decreasing modulus : 
| Ao | = 1 > | Ai| > | A2 1 > • • The Ritz eigenvalues are numeri- 
cally correct (with an error ~ 10 -14 ) for more then 1000 first 
eigenvalues thus demonstrating the very good convergence of 
the Arnoldi method. Third row : The left panel shows the den- 
sity p(A) of eigenvalues in the complex plane, being normalized 
by f p(A) d 2 X — 1, as a function of the modulus |A| for the val- 
ues M = 100, 140, 200, 280. The peak at |A| = 0.02 is outside 
the plot range and has values p(0.02) = 7.7 (M = 280), 8.3 
(M = 200), 9.0 (M = 140), and 10 (M = 100). The right 
panel shows the density p(A) in the region |A| £ [0.58, 1] for 
M = 280 (full spectrum) and M = 400, 560, 800, 1120, 1600 
(partial spectrum). For 400 < M < 1120 only the largest 3000 
eigenvalues and for M = 1600 the largest 5000 eigenvalues 
were calculated by the Arnoldi method and therefore the cor- 
responding densities deviate from the convergent density curve 
at small values of A. 



us to compute a considerable number of eigenvalues (with 
largest modulus) and the associated eigenvectors of S for 



the values M = 400, 560, 800, 1120, 1600 corresponding 
to the matrix dimension of the UPFO Nd = 33107, 63566, 
127282, 245968, 494964 (for the map Q at K g ) which are 
absolutely inaccessible by a full matrix diagonalisation. 
For the case with strong chaos at K = 7 or the separa- 
trix map the matrix dimension is even close to Nd « 10 6 
for M = 1600. In order to provide for a self-contained 
presentation, we give a short description of this method 
here. 

The main idea of the Arnoldi method is to construct 
a subspace of "modest" , but not too small, dimension ua 
(in the following called the Arnoldi- dimension) generated 
by the vectors f , Sfo, • •• # nA_1 £o (called Krylov 
space) where £o is some normalized initial vector and to 
diagonalize the projection of S onto this subspace. The 
resulting eigenvalues are called the Ritz eigenvalues which 
represent often very accurate approximations of the exact 
eigenvalues of S, at least for a considerable fraction of the 
Ritz eigenvalues with largest modulus. 

To do this more explicitly, we first construct recur- 
sively an orthonormal set (of nA+1 vectors) £o> £i > • • • ? £n A • 
For k = 0, 1, . . . , ua — 1 we define the vector Vk+i as the 
Gram- Schmidt orthogonalized (but not yet normalized) 
vector of 5 with respect to £o> • • • ,€k an d store the ma- 
trix elements hj^ — <£j \S \ (,k> for j = 0, . . . , k which 
were used during the orthogonalization scheme. Further- 
more we define the matrix element hk+i,k = 11 ^fe+i || and 
normalize Vk+i by = Vk+i/hk+i,k- Then the product 
S £k can be expressed in terms of the orthonormal vectors 
Zi by: 

fc + 1 

and therefore the matrix hj^ is the representation matrix 
of S in the Krylov space. This expansion is called in the 
mathematical literature [36] Arnoldi- decomposition when 
written in matrix form and it is actually an exact identity. 
However, it is not closed since S £k requires a contribution 
of £k+i unless hk+i,k — for some value of k in which 
case we would have obtained an exact ^-invariant sub- 
space and the diagonalization of the Arnoldi matrix hj^ 
would provide a subset of exact eigenvalues of S (those 
with eigenvectors in the S'-invariant subspace). An inter- 
esting situation appears if due to numerical rounding er- 
rors hk+i,k is ver y small and not exactly zero. Then the 
method automatically generates, with the help of round- 
ing errors, a new "pseudo-random" start vector and ex- 
plores a new subspace orthogonal to the first S'-invariant 
subspace which is actually useful to obtain further eigen- 
values. 

However, when diagonalizing the UPFO S for a chaotic 
map with large dimension this situation, which may be 
quite important in certain other cases, does not happen 
and hk+i,k is always different from zero (actually hk+i,k is 
quite comparable in size to the modulus of eigenvalue A&). 
Therefore we have to cut the above iteration at some max- 
imal value of k. In order to calculate the Arnoldi matrix of 
dimension ua one must actually be careful to determine 
nA + 1 vectors, otherwise one would miss the last column of 
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the matrix h. We also note that the Arnoldi matrix hj^ is 
of Hessenberg form (hj^ = if j > k + 1) which simplifies 
the numerical diagonalization since one can directly call 
the subroutine for the Qi^-diagonalization and omit the 
first, quite expensive, step which transforms a full matrix 
to Hessenberg form by Householder transformations. 

We mention as a side remark that for symmetric or 
hermitian matrices S one can show that the matrix hj^ 
is tridiagonal and the orthogonalization needs only to be 
done with respect to the last two vectors resulting in the 
well known Lanczos algorithm. In principal, the use of an 
exact mathematical property, which may be violated due 
to numerical rounding errors, is somewhat tricky and may 
require special treatment in the various variants of the 
Lanczos method. However, the Arnoldi method always re- 
quires orthogonalization with respect to all previous vec- 
tors and does not suffer from this kind of problem but it 
is also more expensive than the Lanczos algorithm. 

The Arnoldi method requires Nd k c double precision 
registers to store the non-zero matrix elements of S, Nd ua 
registers to store the vectors and const, xn^ registers 
to store hj y k (and various copies of h). The computational 
time scales as Nd n c ua for the computation of S with 
Nd n\ for the Gram- Schmidt orthogonalization procedure 
(which is typically dominant) and with const. xn\ for the 
diagonalization of hj^- 

In the practical applications of the Arnoldi method an 
important point concerns the "good" choice of the initial 
vector £o- It is actually a bad idea to chose a vector which 
is close to the eigenvector of maximal eigenvalue (or other 
eigenvalues) because this would suppress contributions of 
other eigenvectors which we want to retain. A much better 
choice is a random initial (normalized) vector. During the 
Arnoldi iteration the method will automatically suppress 
the eigenvector contributions with respect to the smallest 
values | Xj | and retain the contributions of eigenvalues close 
to the unit circle. If the spectrum has some well-defined 
modest gap between Ao = 1 and the other eigenvalues 
the random initial vector is indeed a very good choice 
and we have used this choice for the case of map at 
K = 7 which we discuss in Section 5. However, at critical 
Kg there is no real gap (see for example the upper right 
panel in Fig. [2| and there is also a considerable number of 
eigenvalues close to unit circle. In this case the inherent 
suppression of small eigenvalues by the Arnoldi method 
may not be sufficiently fast, if \Xk\ k is not small for k close 
to the chosen Arnoldi dimension ua- Therefore we have 
chosen here an initial (normalized) vector obtained from 
an initial number of iterations of S applied to a random 
vector : £o S nini - Random with riini. being the number 
of initial ^-iterations which we have chosen to scale with 
M 2 : n ini . = M 2 /200 (except for the case M = 1600 where 
we have chosen Tiini. = 7000). 

As a first illustration, we have applied the Arnoldi 
method with ua = 1500 to the case of M = 280 and 
Nd = 16609 for the Chirikov standard map at K g , for 
which we were still able to diagonalize the full matrix S. 
In the middle right panel of Fig. [2j we show the modu- 
lus of the difference (in the complex plane) of the Ritz 



eigenvalues and the exact eigenvalues as a function of the 
level number j. The first 1000 Ritz eigenvalues, out of 
1500 in total, are numerically correct with a deviation 
~ 10 -14 entirely due to numerical rounding errors. If we 
only require graphical precision (~ 10 -5 ) there are actu- 
ally 1200 Ritz eigenvalues which are still acceptable. This 
can also be seen in the middle left panel of Fig. [2] where 
we compare the spectrum in the complex plane of the 
full matrix S with the partial spectrum obtained by the 
Arnoldi method. This provides a quite impressive confir- 
mation of the accuracy of the Arnoldi method. For larger 
values of M we have done similar verifications, for exam- 
ple by comparing the Ritz eigenvalues for different values 
of ua or for different initial vectors. Choosing typically 
ua = 1500 or ua = 3000 — 5000 for the largest values 
of M = 1120, M = 1600, we always have a consider- 
able number (at least 500 to 1000) of numerically accurate 
eigenvalues. 

Concerning the (right) eigenvectors, we prefer to deter- 
mine them independently by the method of inverse vector 
iteration which provides numerical reliable (real or com- 
plex) eigenvectors with n\ operations per eigenvector due 
to the Hessenberg form of hj^- Suppose that cp is such an 
eigenvector of hj^ with eigenvalue A, 

r iA — i 

A^j = ^2 h i> k ^ k ' ( 4 ) 
k=o 

then we obtain by Eq. ([3]), the corresponding eigenvector 
tp of S directly from: 

n A -l 

^ = Vk£k - (5) 

k=0 

4 Chirikov standard map at K g 

We first apply the Arnoldi method to the Chirikov stan- 
dard map ([lj at the critical value K g = 0.971635406. The 
Arnoldi method allowed us to obtain a considerable num- 
ber of eigenvalues and eigenvectors of the UPFO S for the 
values M = 400, 560, 800, 1120, 1600 corresponding to 
the matrix dimension N d = 33107, 63566, 127282, 245968, 
494964. We choose the Arnoldi dimension ua = 3000 for 
M < 1120 and n A = 5000 for M < 1600 and we also com- 
pute the first 500 eigenvectors for each case. Even though 
we are not able to calculate the full spectrum for these 
cases, the partial densities of the eigenvalues in the com- 
plex plane, for |A| close to 1, are in a good agreement with 
the full densities obtained for M < 280 as can be seen 
in the bottom right panel of Fig. [2] Below we present the 
most important part of obtained data, more details with 
many eigenstates and high resolution figures are available 
at |3Z| . 

In Fig. [3] we show the decay rates 7^ = — 21n(|Aj|) as 
a function of the level number j (with eigenvalues sorted 
by decreasing \Xj\ or increasing jj) for the case M = 800. 
We note that the first 6 eigenvalues follow quite closely 
a quadratic dispersion law jj « 71 j 2 for < j < 5. 
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Fig. 3. (Color online) Decay rates jj = —2 In ( | | ) versus level 
number j (red crosses) for the UPFO eigenvalues Xj of the map 
(l) at K = K g , M = 800 and N d = 127282. The green curve 
corresponds to the quadratic dispersion law jj « 71 j 2 which 
is approximately valid for the diffuson modes with < j < 5. 




These 6 eigenvalues are actually real, positive and close to 
1. Their corresponding eigenvectors (which are also real) 
extend over the full phase space covered by the chaotic 
trajectory used to determine the UPFO. 

In the first row of Fig. [4] we show the density plots of 
the (right) eigenvector tpo(m) in phase space representa- 
tion (the index m gives the discretized phase space po- 
sition of the cell m) for M = 800 and M = 1600. The 
second row shows |^i(ra)| and |^2(m)| for M = 800. In 
agreement with the ergodic theorem, the eigenvector ^0 
represents a nearly uniform density on the chaotic com- 
ponent. The eigenvectors tpj for 1 < j < 5 (and also for 
certain higher values j < 15 if the associated eigenvalue is 
real, positive and close to 1) correspond to some kind of 
"diffuson modes" with a roughly uniform distribution in 
the angle coordinate x and a wave structure with a finite 
number of nodes in the action coordinate y. For example 
\ipi(m) I is maximal at the upper and lower borders of the 
available phase space and ^i(m) changes sign on exactly 
one curve in between. For ^(^Ol there are three maximal 
curves and two node curves with a sign change of ^2(^1) 
and so on for other diffuson modes. Such diffuson modes 
with quadratic spectrum naturally appears as a solution 
of the diffusion equation 



dp = d_ ( dp 
dt dy V v dy 



(6) 



with boundary conditions dp/dy = 0a,ty = and y « 
1 - r g : pj(y) ex cos(7rjy/(l - r g )), jj « 7r 2 D y j 2 /(l - r g ) 2 
(assuming D y to be constant on the interval < y < 
1 - r g ). 

The structure of the eigenvectors for complex eigen- 
values (or real negative eigenvalues close to "—1") is very 
different and corresponds to the "resonance modes" which 
are typically concentrated (or even localized) around one 
(or a chain of few) resonance(s). This can be seen in the 



Fig. 4. (Color online) First row : Density plot of the eigen- 
vector with eigenvalue A = 1 for M = 800 and N d = 127282 
(left panel) and for M = 1600 and N d = 494964 (right panel). 
In last three rows M = 800 and N d = 127282. Second row : 
Density plot of the modulus of the components of the eigen- 
vectors for Ai = 0.99980431 (left panel) and A 2 = 0.99878108 
(right panel). Third row : Density plot of the modulus of 
the components of the eigenvectors for A6 = —0.49699831 + 
i0.86089756 « |A 6 |e i27r/3 (left panel) and A 8 = 0.00024596 + 
i0.99239222 w |A 8 |e i27r/4 (right panel). Fourth row : Density 
plot of the modulus of the components of the eigenvectors for 
A13 = 0.30580631 + % 0.94120900 « |Ai 3 |e 



and A19 — 
panel) . 



-0.71213331 + % 0.67961609 « |Ai 9 | e 



(left panel) 
(right 



i27r(3/8) 



third and fourth rows of Fig. [4] containing the density 
plots of \i(jj(m)\ for j = 6 and j = 8 (third row) and 
for j = 13 and j = 19 (fourth row). The complex phase 
(fj of Xj = \Xj \ e % lfj for such an eigenvalue represents 
quite well the periodicity of a trajectory with a period 
q if (fj w 27r(p/q) is approximated by a rational num- 
ber times 27r. The fraction p/q represents the position of 
the resonance in the rotation number r. In Fig. [4] we can 
identify cp 6 w 2tt(1/3), <p 8 w 2tt(1/4), (p 19 w 2tt(3/8) and 
as a secondary resonance (close to the main resonance at 
p = 0) ^i3«2tt(2/5). 

The density plots of further resonance modes (with 
complex or real negative eigenvalue) are typically simi- 
lar and they approximately repeat the modes associated 
with the largest resonances but with modified phases and 
decay rates. For example the density of the mode A10 = 
-0.99187524 (for M = 800) with the phase cp 10 = 2tt(2/4) 
is very similar to the density of the mode Ag with y? 8 = 
27r(l/4) representing the resonance at 1/4. Another exam- 
ple concerns the resonance modes close to the resonance 
1/3 (see density plot of the mode A6 in Fig. PO). There is 
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a certain number of higher modes with phases that can 
be written in the form 27r(pi/3 +P2/8) with certain small 
integer numbers pi , P2 . 




Fig. 5. (Color online) Density plot of the complex phase of 
the components of the eigenvectors for A6 = —0.49699831 + 
i 0.86089756 « |A 6 | e i2lv/3 (left panel) and A 8 = 0.00024596 + 
% 0.99239222 « |A 8 | e i27v/A (right panel). Deep blue corresponds 
to either empty cells or phase = — 7r, green to phase = and 
red to phase = tt (M — 800 and N d 127282). 



One may also ask the question in how far the complex 
phases of the eigenvector components carry interesting in- 
formation. As a illustration we show in Fig. [5] for two ex- 
amples the density plot of the complex phase of ipj (m) for 
j = 6 (or j = 8) and M = 800. Even though these modes 
are quite well localized (see Fig. [4| close to the resonances 
1/3 (or 1/4) they still extend, with well defined complex 
phases of i/jj(m), to the full accessible phase space de- 
scribed by the UPFO. In the region with very small values 
of \ipj(m)\ the phase dependence is quite complicated and 
one cannot provide a simple physical interpretation. How- 
ever, in the region of maximal \ipj(m)\ close to the classical 
resonance 1/3 (1/4) for j = 6 (j = 8) one can identify a 
simple structure where the phase is roughly constant on 
a boundary layer outside of each stable region associated 
to the resonance but with different values 27r(//3)+const. 
(27r(Z/4)+const.) for each of the three (four) islands char- 
acterized by the number I = 0, 1, 2 (1 = 0, 1, 2, 3). 

Concerning the first non-zero decay rates one impor- 
tant question is the dependencies of jj(M) as a function 
of M and in particular what happens in the limit M — » 00. 
In Fig. 6j we show 71 (M), for the first non-trivial diffuson 
mode (left panel), and jq(M) : for the first resonance mode 
(right panel), versus M in a double logarithmic scale. In 
both cases jj (M) seems to tend to zero for M —> 00 and a 
power law fit for the range 400 < M < 1600 indicates the 
behavior 71 (M) « 2.36 M -1,30 and j 6 (M) « 389 M -1,55 . 
However, the situation for 71 (M) seems more subtle and 
the curvature, when taking into account the range of all 
values 25 < M < 1600, seems to indicate a transition from 
M~ 2 for small M- values to M" 1 for larger values of M. 
Actually the data can also be well described by the fit: 



7 i(M)«/(M) 



D 1 + C/M 
M l + B/M 



(7) 



with D = 0.245, C = 258 and B = 13.1. 

A physical interpretation of this behavior will be dis- 
cussed in Section 7. Here we only note that the vanishing 
limit liniM^oo lj(M) = is coherent with the observation 
that the relaxation to the uniform ergodic eigenvector is 





Fig. 6. (Color online) The left panel shows the decay rate 
71 (M) of the first "excited" diffuson mode (see second row, left 
panel in Fig.|4| of the UPFO as a function of M (red crosses) in 
a double logarithmic scale. The lower (blue) line corresponds to 
the power law fit 2.36 M -1-30 . The upper (green) curve corre- 
sponds to the fit : f(M) = § X^bJm with D = °- 245 ' C = 258 
and B = 13.1. The right panel shows the decay rate 76 (M) 
of the mode associated to the resonance 1/3 (see third row, 
left panel in Fig. 4j as a a function of M (red crosses) in a 
double logarithmic scale. The (green) line corresponds to the 
power law fit 389 M -1 - 55 . Both power fits were obtained for the 
range 400 < M < 1600 while for the fit with f(M) all values 
25 < M < 1600 were used. 



described by a power law decay of the Poncare recurrences 
in time. 

Before we close this section, we come back to the dis- 
cussion of the complex phases of the eigenvalues. As al- 
ready observed above, it seems that for the resonance 
modes the eigenvalue phases (pj of Xj = \\j\exp(i(pj) are 
close to 2tt x p/q, where p/q is some rational number with 
small values of q. In Fig. [7] we analyze the statistical be- 
havior of (pj/(27r) in two ways. Namely, we show (pj versus 
an index j ordered in such a way that cpj < Pj+\. The 
horizontal lines show the rational numbers of the Farey 
sequence of order 8 (i. e. all irreducible rational numbers 
between and 1 with a maximal denominator 8) obtained 
from a continuous fraction approximation of pj/(2n). One 
can see that these rational values correspond to small steps 
indicating a larger than average probability to find a ratio- 
nal number with small denominator. This feature is seen 
even clearer in the top right panel of Fig. [7] where the 
distribution of (pj/(27r) has well defined peaks at the po- 
sitions associated with the Farey sequence. 



5 Chirikov standard map at K = 7 

We now turn to a particular case of strong chaos at K = 7 
which was previously studied in [16] and by Chirikov in 
[38] . According to [16] the statistics of Poincare recur- 
rences on line y = drops exponentially with time. This 
is also a case even if one takes another line for recurrences, 
e.g. y = 1/2. In the latter case the recurrences are mainly 
determined by a trajectory sticking in a vicinity of two 
small stability islands located on a line y = 0. However, 
in this case the border of islands is very sharp and the 
statistics of Poincare recurrence still decays exponentially 
up to rather long times P(t) ~ exp(— t/r) with a decay 
time r = 23.1 |38]. 
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Fig. 7. (Color online) The top left panel shows the or- 
dered complex phases ipj of the complex eigenvalues Xj = 
\\j\exp(iipj) of the UPFO for M = 1600 and the largest 2000 
eigenvalues (obtained by the Arnoldi method with the Arnoldi 
dimension ua = 3000) as a function of the index j such that 
<fj < Pj+i- The horizontal (blue) lines represent the fractional 
values of the Farey sequence of order 8 corresponding to steps 
in the phase number function. The two bottom panels show the 
same data at with higher resolution in a vicinity of ratios 1/3 
(left panel) or 1/4 (right panel) for the phase rotation number. 
The top right panel shows the distribution of the phase rotation 
numbers obtained by a histogram of bin width 1 / 840 and the 
complex phases of the largest 3000 eigenvalues. The positions 
of the local maxima correspond quite well to the fractional 
values of the Farey sequence. Since the complex eigenvalues 
appear in complex conjugate pairs phase numbers (p/2ir > 1/2 
have been mapped to values below 1/2 by ip/2ir — >> 1 — (p/2jr. 



We have determined the UPFO S at K = 7 for the 
same values as previously (25 < M < 1600) using a tra- 
jectory of length 10 11 which is sufficient due to the faster 
relaxation to the ergodic distribution (as compared to the 
case of critical K g where we used a trajectory of length 
10 12 ). Now the matrix size Nd of S is very close to its max- 
imal value M 2 /2 since with the exception of the two small 
stable islands nearly all cells are visited by the trajectory. 
For M = 1600 we have N d = 1279875 . 

We have calculated the full eigenvalue spectrum of S 
by direct diagonalization for M < 140 (corresponding to 
Nd < 9800) and the first ua = 1500 eigenvalues (and the 
first 500 eigenvectors) by the Arnoldi method for M > 200 
(Nd > 20000). As compared to the case of critical K g the 
necessary time and memory resources are increased due to 
larger values of Nd at given M but on the other hand we 
get reliable eigenvalues for smaller numbers of the Arnoldi 
dimension because the modulus of the eigenvalues decay 
much faster with increasing level number. This can be seen 
at the eigenvalue spectrum in the complex plane shown in 
the left panel of Fig. [8] for M = 140. There are only few 
eigenvalues outside the circle of radius 0.5 and we can 
also identify a clear gap between the first two eigenval- 
ues A = 1 and Ai = 0.8963823322 (for M = 140). The 
density of eigenvalues in the complex plane (normalized 




a(l+bM' ! ) c 



y class =0.0866 



-1 -0.5 0.5 1 



0.005 



Fig. 8. (Color online) The left panel shows the eigenvalue 
spectrum in the complex plane of the UPFO (t — 10 11 it- 
erations) for the map at K = 7 for M = 140 and 
Nd — 9800 = 140 2 /2 as red dots. The green curve shows the 
unit-circle |A| = 1. The right panel shows the decay rate 71 (M) 
as a function of M" 1 for 100 < M < 1600. The upper (blue) 
straight line corresponds to the fit 71 (M) = d + eM" 1 for 
M" 1 < 0.004 resulting in d = 0.0994 and e = 20.3 suggesting 
the extrapolation limit limM^oo7i(M) = 0.0994. The lower 
(green) curve corresponds to the fit 71 (M) = a (1 + 6M _1 ) C for 
M" 1 < 0.01 resulting in a = 0.0857, b = 1370 and c = 0.389 
suggesting the extrapolation limit lhriM^oo7i(M) = 0.0857. 
The black dot marks the decay rate j c \ = 0.0866 found di- 
rectly from the the Poincare recurrences in |38| . 



by / p(X)cPX- 
expression: 



1) can be quite well approximated by the 



p(X) « exp (2.55 - 6.0 |A| - 11.4 |A| 2 



|A| < 0.65 



( 8 ) 

and for |A| > 0.65 we have p(X) < 0.001. This expression 
fits the density for all values 25 < M < 140 for which we 
have been able to compute the full eigenvalue spectrum of 
S. 

In view of the exponential distribution of Poincare re- 
currence times as found in [16,38 , a very important ques- 
tion concerns the limit of the first non-zero decay rate 
71 (M) as M — »> 00. For the case of critical K g with a power 
law distribution we found a vanishing limit of 71 (M) but 
here we expect a finite limit. This is indeed the case as 
can be seen in the right panel of Fig. [8] where we show 
7i(M) as a function of M _1 . A simple fit with two pa- 
rameters 7i(M) = d + eM _1 for M" 1 < 0.004 results in 
d = 0.0994 ± 0.0018 and e = 20.3 suggesting the finite ex- 
trapolation limit liniM^oo 7i(^0 = 0.0994. However, as 
can be seen in the figure, the quality of the fit is not very 
good and can be improved by a more suitable three param- 
eter fit: 71 (M) = a(l + 6M- 1 ) c for M _1 < 0.01 resulting 
in a = 0.0857 ±0.0036, b = 1370 and c = 0.389 suggesting 
the extrapolation limit liniM^oo7i(^0 = 0.0857 which 
actually coincides (within the error bound) with the "de- 
cay rate" j c \ = 2/23.1 = 0.0866 found in [38] from the ex- 
ponential tail of the distribution of Poincare recurrences. 

Concerning the eigenvectors of 5, we mention that 
the eigenvectors for the first mode Ao = 1 represents of 
course the uniform ergodic distribution on the (nearly) full 
phase space with exception of the two small stable islands. 
The eigenvector structure is more interesting for the other 
(non-uniform) modes and in Fig. [9] we show for M = 1600 
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6 Separatrix map at A c = 3.1819316 



In this section, we study the UPFO for a different map, 
called the separatrix map [4], denned by : 




Fig. 9. (Color online) Density plot of the modulus of the com- 
ponents of the eigenvectors for Ai = 0.94665516 (top panel) 
and A 2 = -0.49451923 + ^0.80258270 (bottom panel) of the 
UPFO for the map Q at K = 7, M = 1600, N d = 1279875 
and Arnoldi dimension ua — 1500. 




Fig. 10. (Color online) Increased representation of the regions 
of the two stable islands, close to x — 0.33 and y — (left 
panel) or x = 0.67 and y = (right panel), of the eigenvector 
for Ai = 0.94665516 (see top panel in Fig. [5}. 



the density plots of the eigenvectors for the two modes 
Ai = 0.94665516 and A 2 = -0.49451923 + i 0.80258270. 
One can clearly identify the invariant manifolds and a very 
interesting structure around the stable islands. 



In Fig. 10 we furthermore show zoomed density plots of 
the eigenvector for the mode Ai for the two regions close 
to the stable islands at x = 0.33, y — and x = 0.67, 
y = 0. Both islands cover 125 out of 1280000 cells in total 
with a relative phase space volume being approximately 
125/(1600 2 /2) « 9.7 x 10" 5 that, up to statistical fluctu- 
ations, is in agreement with the result 7.8 x 10 -5 of [38] . 



y = y + sin(27nr) 



A 

2^ 



\n(\y\) (modi) . (9) 



This map can be locally approximated by the Chirikov 
standard map by linearizing the logarithm near a certain 
yo that leads after rescaling to the map ([!]) with an effec- 
tive parameter K e s = A/\yo\ [4]. Therefore the separatrix 
map exhibits strong chaos for small values of \y$\ <C A 
while for larger values of \yo\ we have the typical KAM- 
scenario similar to the Chirikov standard map for small 
or modest values of K. 

For the separatrix map the width of the chaotic com- 
ponent \y\ < yb can be estimated from the condition 
K e R & A/yi « 1 that gives y^ « A. It is known that 
the golden curve with the rotation number r = r g = 
(a/5 + l)/2 = 0.618.. is critical at A c = 3.1819316 [16] 
at which there is a quite large chaotic domain confined 
up to values \y\ < 3.84. Therefore we define the M x M 
cells to construct the UPFO for the phase space domain 
< x < 1, —4 < y < 4. As in the case of the Chirikov 
standard map, we use a symmetry: x -+ x + 1/2 (mod 1), 
y -+ — y, to reduce this range to < x < 1, < y < 4 for 
M x M/2 cells. It turns out that for the separatrix map the 
number of cells visited by the trajectory (of length 10 12 ) 
scales as Nd ~ C^M 2 /2 with Cd ~ 0.78 meaning that the 
chaotic component contains about 78% of the total area 
of the domain (e.g. N d = 997045 for M = 1600). 

As in the case of the Chrikiov standard map the ma- 
trix S is very sparse with small numbers n <C Nd of non- 
zero elements per row (or per column). The average of 
these numbers (with respect to all rows or all columns) 
is k c = (k) w 19 (for M = 1600 and with 12 < « c <19 
for the other values of M). However, n has a very large 
distribution p(tz) depending if we consider the number 
of transitions from (or to) a cell which is either in the 
strongly chaotic range for small y or in the range close to 
the critical curve. This distribution has a power law tail 
p(k) ~ 1/ n 2 for the range n c < n < ft m with a very large 
maximal value n rn ^> k c (e.g. K rn — 2123 for M — 1600). 
We also note that the peak position n p of the distribu- 
tion p(k) is considerably smaller than n Cl e.g. n p = 6 for 
M = 1600. The difference between k c and n p is clearly 
due to the long tails of p{n). These features of the matrix 
S are coherent with the effective value i^ e ff. = A/\y\ of 
the chaos parameter which produces large stretching. This 
property does not create any problems for the Arnoldi 
method and gives only a slight increase of the amount 
of required computer resources (in memory and compu- 
tational time) since the dominant contributions to these 
resources for 1 < tla <C Nd come from terms which do not 
contain n c (see Section 3). 

We have been able to calculate the full eigenvalue spec- 
trum of S for the separatrix map for 25 < M < 200 (279 < 
Nd < 16105) and the first 3000 eigenvalues (and the first 
500 eigenvectors) by the Arnoldi method for 280 < M < 
1600 (31273 < N d < 997045). 
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Fig. 11. (Color online) Density plot of the modulus of the 
components of the eigenvectors for Ai = 0.99970603 (first 
row, left panel), A2 = 0.99828500 (first row, right panel), 
A 3 = -0.99816880 « |A 3 | e* 27r(1/2) (second row, left panel), 



I Ais I e 



2tt(5/13) 



(sec- 



A 18 = -0.73824747 - z 0.66068553 
ond row, right panel), A 20 = -0.80147707 + ^0.58216934 « 
|A 20 | e* 27r(2/5) (third row, left panel), and A 26 = -0.89084450 + 
i 0.42827996 « |A 26 | e* 2?r(3/7) (third row, right panel) of the 
UPFO (10 12 iterations) for the separatrix map at critical 
A c = 3.1819316, M = 1600, N d = 997045 and Arnold! di- 
mension tia — 3000. In the first three rows and the fifth row 
the phase space covers the range < x < 1 and < y < 4. 
The fourth row shows zoomed density plots of the eigenvec- 
tors for A20 (left panel) and A26 (right panel) in the phase 
space range 0.45625 < x < 0.83125 and 2.5 < y < 4. The 
fifth row shows two modes in the strongly chaotic region for 
A 77 = -0.49158867 + i 0.85154001 « |A 77 | e* 27r(1/3) (left panel) 
and A79 = 0.98321618 (right panel). 



As previously the mode for Ao is uniformly distributed 
in the available phase space and therefore is not shown as 
a density plot here. In Fig. [TTJ we show the density plot 
of the more interesting modes Xj for j = 1, 2, 3, 18, 20, 
26, 77, 79 (for M = 1600, N d = 997045 and n A = 3000). 
The first two of these modes (first row in Fig. 11) are 
similar to diffuson modes in the Chirikov standard map 



at Kg that is also confirmed by the quadratic dispersion 



of the associated decay rates jj (see left panel of Fig. 12). 
However, the total number of diffuson modes in the list of 
leading eigenvalues is reduced to only three modes (if the 





5 10 15 20 25 30 
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Fig. 12. (Color online) The left panel shows the decay rates 
7j = — 2 In ( I A j I ) versus level number j (red crosses) where Xj 
are the complex eigenvalues of the UPFO of the separatrix map 
for M = 1600 and Nd = 997045. The green curve corresponds 
to the quadratic dispersion law jj w 71 j 2 which is approxi- 
mately valid for the diffuson modes with < j < 2. The right 
panel shows the decay rates jj(M) for j = 1 (red crosses), 
j = 3 (green open squares) and the eigenvector associated to 
the phase 27r(5/13) (blue full circles, see middle right panel in 
Fig. 11 ) of the UPFO of the separatrix map as a function of M 



in a double logarithmic scale. The upper (black) straight line 
corresponds to the fit 7(M) « 3.48 M~ 71 , the middle (cyan) 
line corresponds to the fit 73 (M) « 1.946 M~ 86 and the lower 
(magenta) line corresponds to the fit 71 (M) « 4.177 M~ 1 203 . 
All fits were obtained for the range 400 < M < 1600. 



uniform mode for Ao is also counted as diffuson mode). 
There are however further diffuson modes characterized 
by real positive eigenvalues Xj close to 1 (e.g. for j = 



16, 17). Actually, if we take out in the left panel of 
the mode for j = 3 (which corresponds to a real 
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Fig. [TTj) correspond well to resonant modes 
>7r(l/2^27r(5/13), 2tt(2/5), 2tt(3/7) and can 



4, 5, 
Fig. _ 

negative eigenvalue, see below), the quadratic dispersion 
law extends even up to the first five modes. 

The four modes A 7 for j = 3, 18, 20, 26 (second and 
third rows in 
with phases 27r(l/2) 
be identified with the resonances at 1/2, 8/13, 3/5 and 4/7 
with 2, 13, 5 and 7 stable islands (we remind that phases 
2tt a and 2tt (1 — a) are always equivalent since they be- 
long to the same pair of complex conjugated eigenvalues 
for < a < 1/2). For Ai§ (second row, right panel) this is 
not very clearly visible since the resonance 8/13 is quite 
small in size as compared to the resonance 3/5 which also 
contributes to this mode. However, the eigenvector com- 
ponents are significantly larger in size at the resonance 
8/13 as compared to the resonance 3/5. In the fourth row 
of Fig. [12] we also show zoomed density plots of the modes 
A20 (left panel) and A26 (right panel) (zoom factor 3.2) 
in order to visualize clearly the fine structure of the res- 
onances. For the mode A20 we see (some of) the small 
islands belonging to the resonance 8/13 even though this 
mode is more maximal at the resonance 3/5 (the other 
way round as for Aig). In short these four resonant modes 
show a similar behavior with phases of the form 27r(p/q) 
as for the Chirikov standard map at K g . 

We furthermore note that the significant properties of 
these 6 modes (2 diffuson and 4 resonant modes) are es- 
sentially determined by the phase space region with y > 2 
(KAM region). We have also identified a few number of 
modes which are determined by the strongly chaotic re- 
gion y < 2. In the fifth row of Fig. 11 we show two of these 
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modes for A77 and A79. These two modes are qualitatively 
quite similar to the two modes shown in the last section 
for the Chirikov standard map at strong chaos K = 7 (see 
Fig. [9| . This again confirms the picture that the separa- 
trix map, at one value of the parameter A, covers implicitly 
various regions with different Chirikov chaos parameters 
K G fi = A/\y\. There are also modes which are quite ergodic 
in the chaotic region and those with a resonant structure 
in the KAM region. The higher diffuson modes (those with 
real positive eigenvalues close to 1) have typically a wave 
node structure in the KAM region, which is quite com- 
plicated due to the two big islands for the resonance 1/2, 
and are simply ergodic (or well extended) in the chaotic 
region. High resolution image files for a selected number 
of these and other modes are available at [37] . 

As in the previous sections we have also studied the de- 
pendence of some of the first non-zero decay rates with M 
and their scaling behavior for M — >• 00. As can be seen in 
the right panel of Fig. fl2j these decay rates (correspond- 
ing to the modes for Ai, A3 and Aig shown in Fig. [TTJ) 
can be quite well fitted (for the values 400 < M < 1600) 
by the power law expressions : 71 (M) w 4.177 M -1- 203 , 
73 (M) « 1.946 M" 086 and j(M) « 3.48 M~ - 71 where 
7(M) corresponds to the mode j = 18 for M = 1600 with 
phase 27r(5/13) and localized at the resonance 8/13 (how- 
ever the level number j of this mode changes with M and 
this is not a fit of "713 (M)"). We note that, as for the 
Chirikov standard map at critical K gi these decay rate 
tend to for M -» 00. 



7 Discussion 

The numerical results for the spectrum and eigenvectors 
of the UPFO presented above clearly show that there are 
modes which relaxation rates 7^0 with M — >• 00. For 
the map |l]) at K g we have 71 ~ 1/M 2 for M 2 < C 2 = 
r c « 0.66 x 10 5 and 71 - 0.2/M for M 2 > C 2 (see Fig.^. 
We interpret this transition in the following way. Accord- 
ing to the results obtained in [18 the average exit time 
r n from an unstable fixed point of the Fibonacci approxi- 
mant r n = p n /q n of the golden rotation number scales as 
r n « r g q n with r g = 2.11 x 10 5 . Thus even for a moderate 
value of qs = 3 we have a very large exit time T3 ~ 6 x 10 5 
which is much larger than I/71 for any M reached numer- 
ically. The Ulam method creates effective noise amplitude 
±1/(2M) in x, y that generates a diffusion with the rate 
Djj rsj 1/(12M 2 ). Due to this noise a trajectory crosses the 
whole interval < y < 0.38 up to the golden curve on a 
time scale tjj ~ 0. 38 2 /Djj ~ 1.73M 2 which is smaller than 
r 3 for M < 600 ~ C. Thus for M < C the smallest re- 
laxation modes have a diffuson type with 71 ~ 1/M 2 . For 
M>Cwe should have T3 <C tjj and the dominance of the 
diffuson modes at low 7 should disappear. There is such 
an indication in Fig. [6] where the crossing between 71 and 
76 should appear at rather large M values. But at the val- 
ues M = 1600 reached in our numerics we only start to see 
an intermediate behavior with 71 ~ 1/M. Thus we think 
that the diffuson modes will disappear at values M > 10 4 
which are unfortunately are out of reach of our numerical 



data. Other modes like 76 correspond to sticking of tra- 
jectories in a vicinity of stability islands. However, it is 
most probable that the lowest values of 7 for such modes 
are also affected by the noise of Ulam method for similar 
reasons as for the diffuson modes discussed above (but on 
a smaller scale around main sticking islands). 

A similar situation appears also for the separatrix map 
where the average exit time T2 ~ 460 from a vicinity of an 
unstable fixed point of the resonance q = 2 is also rather 
large (we determined this time in a similar way as in [18] 
This time is significantly smaller than T3 of the map (ll 
due to strong chaos at \y\ < 2. Due to that we see no 1/M 
behavior for 71 and observe only an intermediate behavior 
1/M which should disappear at larger values of M. The 
modes localized around resonant islands are characterized 
by a decay of their corresponding lowest 7 ~ 1/M 0,8 (see 
Fig. 12 right panel). This dependence on M clearly shows 



that these modes are also affected by the noise of the Ulam 
method. 














jo *fj^ 1-4995 







Fig. 13. (Color online) Rescaled level number j /Nd versus the 
decay rate 7^, in a double logarithmic scale, for the Chirikov 
standard map at K g (left panel) and the separatrix map (right 
panel). Red data points correspond to M = 1600, green to 
M — 800, blue to M = 400 and magenta to M = 200 (from 
bottom to top at jj = 0.2). The black straight line corresponds 
to the power law fits j/N d « 0.052745 7 15203 (left panel) and 
j/N d & 0.014174 7 14995 (right panel) using the data for M = 
1600 in the range 0.04 < 7 < 0.3. The statistical error bound 
of the exponents obtained from the fits is close to 0.1% in both 
cases. Here we used ua = 3000 for the Arnoldi method at all 
M, except the map Q at M = 1600 with n A 5000. 



In fact our aim is to recover the properties of the con- 
tinuous Perron-Frobenius operator using the UPFO as a 
convergent approximant. Our results presented above in 
Fig. [2] clearly confirm this convergence at values |A| < 0.9. 
In Fig.[l3]we demonstrate this convergence even at smaller 
values of 7. Indeed, the data of this figure show that the 
integrated density pu(jj) — j/Nd, which gives the rela- 
tive number of states within the interval [0, 7^], is well de- 
scribed by the dependence ^(7) = Ajjj^ 3 (we remind that 
we order 7^+1 > 7j). The prefactor A^ varies by a factor 
2 when the matrix size ex M 2 is changed by a factor 
64 (when changing M from 200 to 1600). We attribute 
this to the fact that there is a small decrease of effective 
measure near the islands with the increase of number of 
cells. However, this growth is saturated at large M and we 
can consider that Ajj — >• const at M — >> 00. While a small 



variation of A^ with M is visible in Fig. 13 the exponent f3 



12 



K.M.Frahm and D.L.Shepelyansky: Ulam method for the Chirikov standard map 



remains independent of Nd within few percents accuracy. 
For the largest value of M = 1600 we obtain f3 = 1.520 
for the Chirikov standard map at K g and (3 = 1.499 for 
the separatrix map (with a statistical error of 0.1% by a 
fit in the range 0.04 < 7 < 0.3). Thus our results show the 
existence of universal dependence Pe{i) ^ 7 1 ' 5 indepen- 
dent of M. This dependence works down to smaller and 
smaller values of 7 when the size M increases (the lowest 
values of 7^ depend on M as we discussed above). 

Thus our results obtained by the generalized Ulam 
method show that the integrated spectral density decays 
algebraically at small 7: 



pz(l) ~ 7^ » /3 w 1.5 . 



(10) 



This behavior leads to an algebraic decay of Poincare re- 
currences P{t) cx 1/tP. Indeed, the probability to stay in 
a given domain e.g. < y < 1/4 can be estimated as 
P(t) ~ J 1 (^(7)/d7)exp(-7t)d7 - 1/^. The case of 
P = 1/2 corresponds to a diffusion on an interval where 
the diffusion equation ^ gives ~ ir 2 D y j 2 /(l — r g ) 2 . In 
this case j oc ^fy] and we have P(t) ~ l/Vi as discussed 
in p!4| [T6] (for t < I/71). For ^ = 1.5 we have the decay 
P(t) cx 1/t 1,5 in agreement with the data for the Poincare 
recurrences found for these two maps (see pT8| [T9j ) . 

The above arguments give an interesting simple rela- 
tion between the exponent of Poincare recurrences and 
the exponent of the spectral density decay. Of course, our 



numerical data for the UPFO spectrum in Fig. 13 have cer- 
tain numerical restrictions showing the algebraic behavior 
in a moderate range 0.03 < 7 < 0.3. Also it is known 
that the power law decay of P(t) has certain oscillations 
of the exponent /3. Thus further studies of the relations 
between the Poincare recurrences and the spectrum given 
by the generalized Ulam method are highly desirable. At 
the moment, on the basis of our data we make a conjec- 
ture that the both exponents are the same. These points 
will be addressed in more detail elsewhere [ 39] . 

In conclusion, our results show that the generalized 
Ulam method applied to symplectic maps with divided 
phase space converges to the Perron-Frobenius operator 
of the continuous map on a chaotic component. The spec- 
trum of this operator has a power law spectral density of 
states (10) for modes with relaxation rates 7^0. The 



exponent of this power law is in agreement with the ex- 
ponent of Poincare recurrences decay established for such 
maps, even if the range of algebraic decay of the spectral 
density is rather moderate compared to the range reached 
for the algebraic decay of the Poincare recurrences. More 
direct relations between the UPFO and the Poincare re- 
currences require further investigations which are in our 
future plans [39 . 
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